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ABSTRACT 

A new ART code based on the iterative refinement method of least 
squares solution for tomographic reconstruction is presented. Accuracy 
and the convergence of the technique is evaluated through the applica- 
tion of numerically generated interferometric data. It was found that, 
in general, the accuracy of the results were superior to other reported 
techniques. The iterative method unconditionally converged to a solu- 
tion for which the residual was minimum. The effects of increased input 
data error, limited total viewing angle, and reduced number of input 
data were studied. The inversion error was found to be only a function 
of the input data error. The convergence rate, on the other hand, was 
affected by all three parameters. Finally, the technique was applied to 
experimental data and results reported. 
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I. INTRODUCTION 


During recent years, image reconstruction techniques have been 
successfully applied to several fields such as medicine, astronomy, 
electron microscopy, nuclear magnetic resonance, geophysics and optical 
interferometry. An extensive review can be found in Reference 1. 

In 1964, prior to the development of optical holography the 
axially symmetric density field had been reconstructed by using the 
inversion of Abel's integral equation^. Under symmetry conditions the 
Radon transform reduces to the Abel equation. Since 1965, holographic 
interferometry has been applied to aerodynamics, heat transfer and com- 
bustion problems. These developments are discussed in detail by 
Vest (3 ^. 

There are several difficulties associated with the reconstruction 
of phase objects: 

(a) No sharp boundaries can be defined for phase objects. The 
definition of the reconstruction region is based on 

a priori information about the flow field. In general, the 
contribution of the phase object outside the reconstruction 
zone is assumed to be negligible. 

(b) Relatively large errors are present in the interferometry 
data. The optical path differences are recorded as inter- 
ference fringes. The accuracy of the fringe data ordinar- 
ily may not be better than 1/2 of the fringe spacing. In 
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addition, for an unsymmetric density field, the number of 
fringes for each projection varies for different view 
angles resulting in a nonuniform data set. The ideal pro-- 
jection data is a set of non-overlapping, equally spaced, 
parallel rays covering the whole reconstruction region. 

(c) Only the relative phase shift of the projection can be 

measured by the interferogram data. There is, therefore, a 
possibility of incorrect identification of fringe numbers. 
This results in a set of "slightly” inconsistent data. 

The mathematical reconstruction method employed in such cases 
should not be overly sensitive to the "noisy" or inconsistent data. 

Although the Radon transformation 

+» 

£(p,0) * // n(x,y)6[p -xcos(0)- ysin(8)]dxdy (1) 

—CO 

provides a rigorous solution to the problem of reconstruction from 
projections, the solution is uniquely determined only by an Infinite set 
of perfect projections. Here, l is the optical path length, n is the 
refractive index, 6 is the Dirac delta function, 6 is the projection 
angle, p is the coordinate along the projection plane, and x and y are 
the coordinates describing the reconstruction region. In practical 
problems, the discrete nature of projection data and the unavoidable 
measurement errors may result in the failure of reconstruction. 
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In this paper, the iterative refinement method of least squares 
solution for tomographic interferometry is discussed. A complete set of 
computer codes for density reconstruction from holographic interfero- 
grams was made. The effects of projection number, limited viewing angle 
and the measurement error on the reconstructed image are determined 
through numerical experimentation. The flow field around the tip region 
of a revolving helicopter rotor blade was reconstructed from 40 inter- 
ferograms. The results of the reconstructed density field are pre- 
sented. 


II. ITERATIVE REFINEMENT METHOD OF LEAST SQUARE SOLUTION 

In principle, the ART (algebraic reconstruction technique) algo- 
rithms are the schemes for solving a large system of linear equations. 
The reconstruction region is divided into a square grid (M =* m x m cell) 
and the refractive index within each cell or pixel, n^ j , is assumed con- 
stant. The Radon transformation then reduces to a set of discrete 
linear equations 

m m 

l l w (p,0) n - L(p,0) (2) 

i-1 j=l 1J 1J 

where W^j(p,9) are weight factors determined from geometric relations. 

It should be noted that only those factors associated with pixels 
through which the projection ray passes are nonzero. For k different 
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projection angles, each having n nonoverlapping and equally spaced 
parallel rays, the number of equations N, is given by N * k x n. 

Equation (2) can be written in the form 

C X + L * 0 (3) 

where, L and X are N and M dimensional vectors respectively, and C is 
the coefficient matrix with N X M elements. Equation (3) can be trans- 
formed into a symmetric, positive definite, NORMAL EQUATION given by 

C T L + C T C X - 0 (4) 

In practice, due to measurement errors and other inconsistencies. 
Equation (3) becomes 

R - C X + L (5) 

where R is the residual vector. An approximate solution is sought for 

which the Euclidean norm of the residual vector HRl is minimum. The 

Euclidean norm is defined as IRI “ vT R. • R . • The solution of the 

i 1 1 

normal equation gives the standard least square solution of (3). 

Procedures for obtaining solutions to a large system of linear 
equations can be found in a number of different scientific areas. Iter- 
ative techniques have been widely used. It begins with an initial esti- 
mate and then it repeatedly modifies the estimate until some threshold 
condition is satisfied. There are different ways to modify the 
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estimate. The effectiveness of a method depends on the matrix 
character. Considering the features of matrix C (large and sparse) and 
relatively large measurement errors, the method of CONJUGATE GRADIANTS 
is found to be applicable for image reconstruction. 

The computer procedures of conjugate gradiant method are as 


follows: 

Initial vector 

k - 0: X(0) - 0; R(0) - L (6) 

For iterative number k * 1,2... 

r(k-l) = C T R(k-l) (7) 

s(k) = - r(k-l) for k - 1 

s(k) =* - r(k-l) + s(k-l)lr(k-l)l/lr(k-2)lfor k > 1 (8) 

q(k) = C s(k) (9) 

X(k) = X(k-l) + s ( k ) n r ( k— l)B/ltq(k)ll (10) 

R(k) = R(k-l) + q ( k ) II r ( k— 1 ) B /lq(k)J (11) 


The matrix C has only five percent nonzero elements. It can be 
stored permanently or, in the case of limited storage capacity, may be 
computed according to a simple geometric relation for every ray. The 
vectors that must be stored are X, R, s, and q. Accordingly, a proce- 
dure written in PASCAL running under CPM-86 with MT+86 compiler has been 
written. This procedure runs on IBM PC microcomputer, and is limited to 
20 x 20 mesh size. 
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It can be proven^^ that: The norm of the residual vector R 

decreases with increasing iteration number* Thus, in principle, this 
code is a technique minimizing the residual R to the level determined by 
the inherent errors associated with the computer truncation and 
measurement data. Therefore, the final reconstruction accuracy is 
determined by data noise. Here the errors associated with the 
measurement of the optical path length (fringe lines) are the major 
attributor to the data noise. 

To evaluate the characteristics of the developed tomography code, 
a numerically generated density field. Figure la, was used as the input 
for the numerical test. It represents the calculated air density dis- 
tribution over the tip region of the rotor blade in a plane above the 
blade The distribution of the refractive index variance, 

x ” (n Q -n), is calculated and shown in Figure lb. The data for the 
integral optical length difference as a function of ray parameters, 9 
and p, was generated from Equation (2). The errors associated with 
optical length measurement, encountered in real tomography, were 
artificially added to the fringe data file. Fringe round-off errors of 
±0.5, ±0.05 and ±0.005 were used here. The reconstruction code was 
parametrically evaluated for different numbers of equations and total 
viewing angles. The numerical test matrix is shown in Table 1. 

Figure 2 shows the reconstructed field for Case 1 after different 
iterations. The overall features are apparent after two iterations and 
further refinement of the solution is obtained with additional iterations. 
The total reconstructed error, R f , is also noted on Figure 2. Here, the 
reconstruction error is defined as 
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( 12 ) 


R r m yj\ E ^ X lj,in “ X ij,reJ 7 \ ^ t x ij,iJ 2 

where and *£j trec are input and reconstructed refractive index 

variances, respectively. 

The dependence of the residual vector norm on the number of iter- 
ations for Cases 1 to 7 of Table 1 is shown on Figure 3. Figures 4, 5 
and 6 demonstrate the effects of fringe number round-off error and view 
angle range on the reconstruction error. From these and other numerical 
experimentations, the following conclusion may be drawn: 

• The convergence rate is increased as the projection ray number, 

N, is Increased. The corresponding Inversion error for a suffi- 
cient number of iterations is, in general, independent of N. 

• The normal of the residual vector, R, in all cases, uncondition- 
ally converges to a limit which is a strong function of the input 
data error. The resulting error for the first few iterations is 
independent of the input data error. However, for more accurate 
input data, additional iterations result in more accurate 
results. This underlines the importance of the accuracy of the 
data set. It also proposes a guideline for determining the maxi- 
mum number of iterations based on the rate of change of the 
residual. 
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Finally, limited look angles, (>30°) do not have a significant 
effect on the convergence rate or the accuracy of the recon- 
structed data. Acceptable results were obtained for a total view 
angle of 30°. Care must be taken to ensure that sufficient 
information is recorded on the interferograms when the view angle 
limitations are imposed. 


III. THE PROCEDURE FOR PROCESSING EXPERIMENTAL DATA 


Figure 7 shows the holographic interferometry arrangement at the 
AVRADCOM Laboratory of NASA Ames for studying the flow around a 
revolving helicopter blade. The coordinate system (x,y) is wing- 
fixed (x-s panwise, y-chordwise). The projection of point "A” on the 
wing defines the origin for coordinate p along any azimuthal angle. The 
photograph of a typical interferogram is shown in Figure 8. 

Forty interferograms were digitized using a Tektronix 4112 dis- 
play terminal and a graphic tablet. The data were reduced and a map of 
the optical path as a function of height above the blade was developed. 
Data at each elevation and for all viewing angles were grouped (Z-files), 
and constituted the input data for the reconstruction code. 

The reconstructed density distribution, (p/p Q ), produced by the 
tomographic reconstruction code is shown in Figures 9 & 10 for the 
heights z = 0.5, 1, 1.5 and 2 inches (chord C - 3 inches; the leading 
edge of the blade Y/C = -0.5; aspect ratio = 13.7). The optimum number 
of iterations was about 3 to 6 for the experimental data. The 
computation time on an IBM PC was approximately 3 minutes per iteration. 
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IV. DISCUSSION 


A new ART code based on the iterative refinement method of least 
square solution is developed for tomographic reconstruction of the 
three-dimensional density field. The features of the technique was 
investigated using numerically generated interferometry data. It is 
shown that, in general, the ART code used here, is capable of recon- 
structing phase objects using limited viewing angle or using data with 
relatively large errors. 

In comparison to other reconstruction techniques, it is believed 
that it is possible to obtain more accurate results with the ART code. 

For similar numerical tests, the reconstructed mean-square error 
for view angles of ±90° and ±40° were reported to be 6.7% and 11%, 
respectively, when the convolution back projections technique was 
used. The corresponding values for R r become 26% and 33%, 
respectively. Similar situations were examined in Cases 3 and 10 with 
R r at only 3% and 10% for 20 iterations. On the other hand, the 
computation time for iterative techniques are generally higher than some 
other techniques. 

An important feature of the technique is that with limited reso- 
lution (say, 20 x 20), it is possible to repeatedly use the same code 
and zoom into the regions surrounding important features of the flow, 
and hence, reconstruct the flow field in greater detail. For this, the 
input data for the subsequent analysis must be corrected. 

Finally, a strong feature of the reconstruction code is its abil- 
ity to accept data within limited view angles. The density field of the 
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numerically generated data was reconstructed for a total view angle of 
30° and with a final error of 10%. This result is extremely important 
when the application of tomography to the flow field within a wind 
tunnel is considered. 
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TABLE 1. THE NUMERICAL TEST MATRIX 


CASE 

FRINGE NUMBER 
ROUND OFF ERROR 

VIEW ANGLE 
(degree) 

NUMBER OF 
VIEWS 

NUMBER OF 
POINTS PER VIEW 

i 

0.005 

0 

• • 

180 

61 

20 

2 

0.005 

0 

• • 

90 

61 

20 

3 

0.005 

0 

• • 

180 

31 

20 

4 

0.050 

0 

• • 

180 

61 

20 

5 

0.050 

0 

• • 

180 

31 

20 

6 

0.500 

0 

• • 

180 

61 

20 

7 

0.500 

0 

• • 

180 

31 

20 

8 

0.005 

-75 

• • 

75 

31 

20 

9 

0.005 

-60 

• • 

60 

31 

20 

10 

0.005 

-45 

• • 

45 

31 

20 

11 

0.005 

-30 

• • 

30 

31 

20 

12 

0.005 

-15 

• • 

15 

31 

20 

13 

0.050 

-75 

• • 

75 

31 

20 

14 

0.050 

-60 

• • 

60 

31 

20 

15 

0.050 

-45 

• • 

45 

31 

20 

16 

0.050 

-30 

• • 

30 

31 

20 

17 

0.050 

-15 

• • 

15 

31 

20 
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2. Reconstructed field of refractive index variance (n Q -n) for Case 1. 
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6. Inversion errors vs. iteration number for different view angles 
(fringe number round off error ■ 0.05 X). 
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10. Density ratio profiles (blade length B * 41.1", chord length 
C = 3"). 


449 


Mi-12/21 


i /uuv 

/-~+Z-~r : L i 

/CZzfrQ W\UU\\nA\K 
/-Z£?-?3I l\l\\ \UWA VSSZT'Nw 
V \ \ \ \\ \ \ \ \ 



^TTTSA \W\W 
r _. j^Tt- 


Figure 1. Numerically generated density and refractive index distribution. 
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Reconstructed field of refractive index variance (n -n) for Case 1. 
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ITERATIVE NUMBER 

Figure 3. Residual vector norm vs. iterative number for Cases 1-7. 
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Figure 4. Inversion errors vs. iterative number for Cases 1-7. 
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Figure 5. Inversion errors vs. iteration number for different view angles 
(fringe number round-off error “ 0.005 A). 
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Figure 7. Rotating blade tip region and the coordinates. 
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Interferogram of flow over a rotating blade, 0 == 176, 
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Figure 9. Reconstructed density fields for different heights above 
blade chord line. 
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